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Abstract. We present a numerical technique for solving evolution equations, as 
the wave equation, in the description of rotating astrophysical compact objects 
in comoving coordinates, which avoids the problems associated with the light 
cylinder. The technique implements a fast spectral matching between two domains 
in relative rotation: an inner spherical domain, comoving with the sources and 
lying strictly inside the light cylinder, and an outer inertial spherical shell. Even 
though the emphasis is placed on spectral techniques, the matching is independent 
of the specific manner in which equations are solved inside each domain, and can 
be adapted to different schemes. We illustrate the strategy with some simple but 
representative examples. 
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1. Introduction 

A general motivation for undertaking this project is suggested by the numerical study 
of astrophysical compact objects, e.g. neutron stars and black holes, individually 
rotating or in bound binary systems, where the dynamics of the interaction field is 
described by a system of equations which is (at least, in part) hyperbolic. An example 
of this is provided by Einstein equations in General Relativity or the Maxwell equations 
in modeling the pulsar magnetosphere. The main objective of this paper is to present 
a general numerical technique for solving the dynamics of such systems, with the 
example of the scalar wave equation as the equation governing the evolution of the 
system. 

The study of the general motion of compact finite bodies, say black holes (BH) 
or neutron stars (NS) for concreteness, in a fixed numerical grid represents a serious 
challenge from a numerical point of view. Examples of this Eulerian approach to 
the motion problem can be found in Refs. p] E]- Recent examples of moving BH 
in the context of binary systems have been successfully developed an implemented 
both using excision techniques [3 [H [5] and moving punctures [6j [3 [8]. However, 
we would find unfortunate if these extraordinary achievements actually shadowed the 
development of parallel approaches which employ independent analytical formulations 
and/or numerical methods. This is the case of constrained evolution formalisms (e.g. 
see Ref. [5]) and of the use of spectral methods (already initiated in the binary 
BH case in Ref. [TO])- More concretely, the use of spectral methods seem to favor an 
excision approach for BH evolution where the coordinate location of the excised surface 
do not change significantly in time [10 . On the other hand constrained evolutions 
schemes, in combination with a comoving coordinate systems, can be employed to 
endow the excised surface with a geometrical character by enforcing appropriate inner 
boundary conditions in some of the elliptic (constraint) equations, something that 
can be relevant both for the numerical stability of the code and the study of some 
quasi- local geometrical properties of BH space-times [11] . These elements indicate the 
interest in pursuing the study of corotating coordinate systems. 

Regarding the presence of matter, namely the NS case, experience in our group 
shows the convenience of using a Lagrangian scheme to describe the surface of the 
NS and a Eulerian one for the interior. Details of this technique can be found 
in Refs. [12] and in |13j where a steady state approximation was used. Building 
accurate non-stationary numerical models of magnetized rotating neutron stars and 
of the surrounding magnetosphere might thus also be easier with such codes, able to 
naturally follow the surface of the star. Therefore, an interesting alternative strategy 
to the Eulerian approach consists in reducing the motion of the body in the grid 
by choosing a comoving coordinate system. However, the use of a single rotating 
coordinate system faces the problems coming from the superluminal motion of the 
grid at distances to the rotation center larger than the light cylinder radius (see next 
section for the definition of the latter). Part of these problems, related to the tensorial 
nature of the rotating fields, can be handled by using appropriate techniques (e.g. the 
"dual-coordinate frames" introduced in Ref. [TO]). 

The aim of this paper is to propose a complementary approach that bypasses the 
light cylinder issue in the case of a scalar field wave equation, using an intermediate 
strategy in which space is split into two domains: first, an internal spherical one 
of radius i?» containing the sources and set inside the light cylinder, and then an 
external domain with inner boundary given by the sphere r = and whose outer 
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radius is larger than the considered wavelength in order to impose approximate (or 
exact) outgoing wave conditions. The exterior domain is described in terms of non- 
rotating coordinates and simultaneously, a rig idlj§ rotating coordinate system is set in 
the interior domain whose angular velocity is fixed by a physical scale in the problem 
(e.g. the orbital frequency in a binary problem). Each domain could be further 
decomposed into smaller ones, but this paper focuses on the matching conditions and 
not in the specific manner in which the equations are solved inside each domain. The 
main problem we address here is then the matching technique between domains in 
relative rotation in a multi-domain spectral scheme. In particular, if we consider the 
spherical boundary between the inner and outer domains, and denote by Ng and N v 
the number of sampling points in the 9 and ip spherical-coordinate directions, then the 
technique for matching the solutions at this intermediate boundary is said to be fast in 
the following sense. The most straightforward approach is just to write the matching 
system, using a spectral summation, at every point of each grid, which requires NgN^ 
operations. Our approach here can be performed with only Ng ■ N v instead, with a 
so-called stroboscopic technique based on spectral methods. 

The paper is organized as follows: Section [2] describes the simple model we want 
to study and specifies the light cylinder problem, Section [3] gives the numerical details 
of our matching technique and some outlines of the implementation, Section [4] shows 
the results of two test-problems, and finally Section [5] contains concluding remarks. 
We emphasize that in this work we address the description of a technique enabling 
us to study the dynamics of a system with a rotating source. Even though the main 
problem we have in mind is the study of the evolution problem for some given initial 
data, we can also employ these resolution techniques in order to construct the initial 
data themselves by assuming additional hypotheses and/or approximations for some 
given terms in the equations. In particular, in the description of the present technique 
we do not make any assumption about the existence of approximate Killing vectors, 
something that would have to be incorporated in the equations themselves in the 
construction of initial data of binaries in quasi-circular motion (see Refs. [HI IT6 l I17 j ) 
for different approaches, and below). 

2. Light cylinder and simplified mathematical models 

The resolution of the full Einstein Equations (EE) can be reduced to solve a couple 
of wave equations for two independent scalar potentials if one uses a constrained 
scheme, where the ten EE have been decomposed as four elliptic constraint equations, 
four gauge degrees of freedom and two scalar wave-like evolution equations [9]. In 
this section we consider the linearized version of the actual problem for one of the 
potentials. This is quite general since in the overall problem for the solution of the 
EE, we shall perform exactly the same matching for every wave-like equation. 

In particular, consider a non-rotating coordinate system (t, r, 9, (p) and a field $ 
satisfying 



where c is the propagation speed of the field and n denotes the matter mass density 
(renormalized with a AttG factor). We assume that matter is contained in a compact 

% The notion of "rigidity" employed here, is denned in terms of a fiducial time-independent flat 
3-metric fij, as introduced in the constrained evolution formalism presented in Ref. [9], 




(1) 
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region bounded by characteristic radius Defining a spherical coordinate system 
rotating rigidly with constant angular velocity (t' =t,r' = r,6' = 6,(p' = (p + Qt), 
and writing explicitly the Laplacian, the wave equation becomes 



We note that in our formulation of the problem, there appears a second characteristic 



cylinder of radius Rl = with respect to the rotation axis, where the vector d t > 
(corresponding to d t — ficL in the inertial system) becomes null. The total operator 
in Eq. §Z§ is hyperbolic everywhere. However we note that, even though the spatial 
part of the operator is negative definite in the region inside the light cylinder, r < Rl, 
this feature is lost for r > Rl- Motivated by this remark (see discussion below), we 
will always assume the hypothesis that i?* < Rl- This assumption will be eventually 
justified on physical grounds. 

As mentioned above, in the present work we study the dynamical evolution of 
a system with a given matter source n(t, r, 8, <p). Then, our strategy to describe the 
dynamics consists in solving the equations in different coordinate systems depending 
on the domain. In concrete terms: a) we choose rotating coordinates in the domain 
containing matter, i.e. we solve Eq. Q for r < and b) we employ an inertial 
coordinate system in the exterior region, i.e for r > R* we solve Eq. ([1]). 

Even though we are focused on the general evolution problem, it is convenient 
to highlight that the existence of the light cylinder is an important issue in the 
construction of quasi-stationary rotating configurations of binary systems. The 
introduction of a helical Killing vector has been used in the literature (cf. p~5| . [TBI ITT] ) 
in order to model slow-motion adiabatic configurations of binary NS and BH systems. 
In particular, regarding Eq. ^ this imposition of the vector dt' as a helical symmetry 
of the solution, implies the vanishing of the first two terms in the left-hand-side of the 
equation. A discussion of the numerical technical issues of the resulting mixed-type 
partial differential equation, that is elliptic inside the light cylinder and hyperbolic 
outside, can be found in the series of works on quasi-stationary binary inspiral 
[P5I [T51 HZUl l2"Tj . as well as the so-called "periodic standing- wave" approximation 
p?!Zl EMI HZ5] , See also Ref. [55] for a study of its mathematical well-posedness. 
In this context of potential numerical problems associated with the type of the spatial 
part of the differential operator, the advantage of our approach is that the light cylinder 
problem is completely avoided, since i?* < Rl- Indeed, in the interior domain the 
differential operator presents a standard "negative" spatial part and there are no 
problems in solving numerically Eq. @, whereas no light cylinder problem shows up 
in Eq. |T]) in the exterior domain. 

We conclude the section by emphasizing that with this point of view, the light 
cylinder problem is replaced by a matching problem, which is much easier to treat 
numerically. We therefore state the main problem to be addressed in this approach, 
and which defines the goal of the paper: the implementation of the matching between 
the solutions in both domains. Next section presents a fastEl matching numerical 

+ By fast we mean that the number of operations is not larger than N ■ logTV for each dimension, 




- df., - —d r 



(2) 





This defines the light 
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algorithm. 

3. Numerical strategy 

We introduce the following dimensionless variables £ and r in the inner and outer 
domains: 

r' = l - (for r' < R,), and 

£ = ^, r = Jj (for r > J?*), (3) 

where P = 27r/f2 is the rotation period. Making use of these coordinates and denoting 
the field <& by < i> < in the interior domain (0 < £ < 1), and by in the exterior 
domain (£ > 1), we write Eqs. |T]) and ([2|) in the following dimensionless form 

(—- C 2 A] $> = , (4) 



\dr 2 



$< = n (5) 



where 



C = 27& (6) 
-ft* 

is the dimensionless "light velocity" , A is the Laplacian in terms of the dimensionless 
radial coordinate £, and n = R 2 C 2 ■ n is the dimensionless matter density (note that 
the dimensionless angular velocity is 2tt). 

In order to simplify the presentation, and since the notation should be self-explanatory 
in the following, we shall drop the "prime" for the equations in the rotating grid 
hereafter. 

3.1. Stroboscopic matching 

We look for a solution of the above Eqs. ^ and ([5]) by expanding the fields $ and n 
in spherical harmonics 

$ = -P™(cos 9) [ai m cos(mip) + b lm sin(m^)] , (7) 

t III 

n = Yj P™ (cos 0) \n\ m cos(mip) + n\ m sm(m<p)] . (8) 

where the coefficients ae m {£,T~) and T ) are the unknown functions, and P™(&) 

are the Legendre associate functions. Similarly, n\ m and n\ m are the components of 
the source n. In the inner domain (£ < 1), we obtain a system of two coupled partial 
differential equations for each value of £ and m (recall that we now use r = r'): 

d 2 c2 ( d 2 , 2 d i{i + 1) 



dr 2 \di 2 e 

d_ 

Ft 



(2^) 2 m 2 



< m + 2(27r)m— 6< m = n c tm , (9) 



where N denotes the number of degrees of freedom in a given dimension. As it has been said, in the 
present case the number of operations is proportional to Ng ■ N v . 
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dr 2 



C 2 



2d_ 

I'd! 



{2Tr) 2 m 2 



£(1+1) 



6< m -2(27r)m— a 



Similarly, for the external domain (£ > 1) we find: 



dr 2 
dr 2 



-C 2 



-C 2 



2d_ 

lai 
2d_ 



e(t + i) 



£(l+l) 



'fni 



= 



= 0. 



(10) 



(11) 



(12) 



By taking into account the expression of the "light velocity" given by Eq. © and 
the inequality m < £, we can see that the second-order differential spatial operator 
appearing in Eqs. J9}, (flOl) is of definite type if Rl > R*- 

In order to solve these equations we make use of a multi-domain spectral technique 
in the radial direction and finite differences in the time variable. In our specific 
implementation, for each couple of values (£, m) and in each domain we construct 
a particular solution as well as the relevant homogeneous solutions; in our simple 
case, regularity of the solution at the origin leaves a single homogeneous solution in 
the interior domain |14j . whereas in the exterior domain there are two homogeneous 
solutions. Solutions to Eqs. ©-(dUJ) in the inner domain and to (TTTjl - lTi"!?]) in the 
outer domain are written as linear combinations of the calculated homogeneous and 
particular solutions, thus involving in total three coefficients to be fixed. The latter 
are determined by: i) imposing an outgoing Sommerfeld wave condition on $ > at the 
outer boundary of the exterior domain (see e.g. Ref. [27]). 

dr d£ £j \^ fD 
where £d is the external radius of the outer domain; ii) enforcing the continuity of the 
fields <I >< and < & > , and Hi) of their radial derivatives $> r< and at the intermediate 
boundary £ = 1. This makes three conditions for three homogeneous solutions and 
the matching system is invertible. 

The following discussions focus on the two last points, i.e. the matching conditions 
at £ = 1, for a given value of the angle 8i, i S {1, Ng}. In our scheme we employ 
the same number of sampling points in the ip direction, N v , in the two grids. Using 
an index j £ {1, ...,N V }, let us denote by 3?J~(t) and $' 3 - (r) the values of the field 
$ < (r, £ = 1,8 = 0i,Lp = ipj) and its radial derivative 3> /<: (t, £ = 1,8 — 8i,ip — ipj) 
on the (^-sampling points in the inner grid, and by <$><"(t) and $'^(7-) the analogous 



1,1 



ipj) and ^(t,^ = 1,1 



, (p = ipj) in the outer 



values of $ > (r, £ 
domain. 

Then, at the initial time r = r Q , we choose the sampling points of the two grids 
to coincide. In this situation we use strong matching technique^, that are simply 
expressed as: 

$<(r )=d>>(r ) , $'<(ro) = $f(ro) 



^(to) 



(to] 



(to) 



(to) = 



* Alternatively, variational (weak) matching methods could be used, 
introduced if this is suggested by the complexity of the system. 



(14) 
(to) 

Penalty terms can also be 
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Figure 1. Choice of evolution step St, adjusted to the sampling interval 
Sip, so as to guarantee the (stroboscopic) coincidence of the ip— sampling points 
corresponding to the interior and the exterior grids along the whole time evolution. 



This makes 2N V conditions for each value 9i, which is the right number of matching 
relations. Regarding the choice of the evolution time step St, two situations may 
occur: if one uses an explicit time-evolution scheme, it must satisfy the stability 
Courant conditions; on the other hand, if the time-integration scheme is implicit 
(more precisely, A-stable), then the time-step is not limited. Let us first suppose 
that this is the case and let us choose St = where Sep — jf- is the sampling 
(^—angular interval. With this choice of St the sampling points at r = r Q + k ■ St, with 
k € {0, 1, 2, ...}, do coincide again at any later time (see Fig. [1]). The matching for a 
step k is simply given by: 



<J>< = <J>> <b'< = <J>'> 

<J>< — <J>> 8>'< — <$>'> 

2 * (2+fc)mod(A^) ' Y 2 ^(2+fc)mod(AT<p) 

3>< — S>> <b'< — (£>'> 



(15) 



We now consider the more general case in which St is bound to be smaller than 
a given critical value, given for instance by the Courant condition. Our strategy 
consists in choosing a St, satisfying St — -^^Sip where if is a sufficiently large 
integer number. In this case, the sampling points only coincide after K time steps, 
and therefore we need to interpolate at supplementary points in the ip coordinate (see 
Figure [5]). At this stage, spectral techniques provide an straightforward manner to 
proceed by performing an "oversampling" of the spectral description. That is, given 
the original N v coefficients codifying the behavior of the function in the ip direction 
and corresponding to the N v sampling points [e.g. {an, an, ■■■,aw v ) for a given £], 
one increases the number of sampling points in ip by a factor K, in such a manner that 
the new coefficients, counterpart of the new points, are set to zero: the information is 
encoded in the same number of non- vanishing coefficients but in a space with a bigger 
number of degrees of freedom: 

n v times 

(aii,ai2,--,at N 0, ...,0 , 0, 0). (16) 

V v ' 

(K-l)-N v times 

The values at the new (interpolation) points can be obtained by performing the 
inverse Fourier transform on the larger set of coefficients. Moreover, a Fast Fourier 
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Transform algorithm could be used only involving a number of operations proportional 
to (N v ■ K)\og(N v ■ K) for each £. Once the matching is done, in an analogous manner 
to the one described in conditions (fTS)) . only N v coefficients are kept in order to 
codify the function $ at each domain. Proceeding in this way, the functions $<, 
and < i> > , can be interpolated at the Gauss-Radau sampling points of the 
Fourier expansion. These Gauss-Radau sampling points are "smart" sampling points 
for which the following identity holds 

[■2k n-i 

/ cos(fc0) cos(ni(j))d(j) — cos(2irkj /N v ) cos(2irmj / N v ) 

fc, to, j being integer positive numbers and k, m < N v /2 (analogous relations for 
other trigonometric functions). If the solution of the above system of equations does 
not contain spatial frequencies (in the Fourier expansion) larger than N^/2, then a 
solution obtained with a number of degrees of freedom higher than N v will give the 
same result (within round-off errors); therefore the interpolation is exact. 

Even though the previous discussion illustrates the underlying philosophy, the 
method presented above is not the most efficient one. As a matter of fact, there is a 
manner to by-pass the inverse Fourier transform step, gaining a factor logiV in each 
angular direction. This is accomplished by making the matching in the coefficient 
space rather than in the configuration values of the functions $. For this, we must 
expand the inner and outer solutions in the same ^-coordinate expansion bases. This 
can be achieved straightforwardly, since the respective bases [cos(mip') and sin(mip') 
in the inner domain, and cos(mip) and sin(mip) in the outer domain], are related 
by a simple rotation: Sip = ip' — ip — flSt — 2hSt. Therefore, given the function 
F(ip') = $ < (r, £ = 1, 9, if') in the inner domain with 

F ( ( f')= ( a ™ cos ( mi P') + b 'm sm(m(p')) , (17) 

m=0 

the standard trigonometric relations provide the coefficients a m and b rn in the ip bases 
employed in the outer domain, and necessary for performing the matching: 

a m = a' m cos(to Sip) + b' m s'm(m Sip) 

b m = — a m sin(TO Sip) + b' m cos(to Sip) (18) 

If Sip = 2tt/N v> this is equivalent to sample the function F(ip) at the Gauss-Radau 
points and to come back in the Fourier space. In other words, the matching is done 
without leaving the Fourier space, but at the cost of (again) oversampling points in 
the (/^-direction by a factor K. The number of operation is proportional to N§ ■ N v . A 
bonus of this strategy, that we refer to as stroboscopic, is the following: the technique 
described above prevents the appearance of spurious time frequency terms generated 
by the beating between the rotation frequency and the time sampling frequency. This 
terms could spoil the study of the numerically constructed signal, in e.g. an a posteriori 
time Fourier transform analysis of the emitted radiation. 

3.2. Implementation 

Equations (|5|)- (fTUl) and (fTT|) - (fT^|) are solved by using the spectral methods developed 
in our group [25] . 
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Figure 2. Choice of evolution step St as an integer fraction (in this case K = 3) 
of the ip— sampling interval Sip. In this case, interpolation of the points is needed 
for enforcing the matching. This is achieved in a fast manner by employing a 
spectral oversampling technique based on the Gauss-Radau points. 



For the time evolution in the non-rotating domain [i.e. Eqs. (fTT| and [)12p] we 
use a second-order implicit Crank-Nicolson scheme. In addition, Sommerfeld outgoing 
boundary conditions are imposed at the outer boundary. 

Regarding the rotating domain, which contains the coordinate singularity at 
£ = 0, we perform an expansion on a Galerkin basis satisfying the regularity conditions 
|14j . Moreover, the wave operator is also treated implicitly with a second-order 
Crank-Nicolson scheme and the coupling terms containing first-order time derivatives 
are treated explicitly by making a second-order extrapolation from known values at 
previous time steps. An implicit treatment of those terms could be more appropriate 
but not necessary. 



4. Tests 



We perform two tests of the scheme previously described. 



L = ^ + ~--^ + (-) . (19) 



4-1- Comparison with an analytical solution 

The first example deals with the comparison between the numerical and the analytical 
solutions, in the restricted case of a helically symmetric operator obtained from Eq. 
© or Eq. ([10} by setting to zero the time derivatives. The resulting differential 
operator, L, has the following form 

d 2 2d 1(1 + 1) t (2m> 2 

de + ldt e 

A straightforward manner to proceed consists in obtaining the source n by the 
application of the operator L to a potential $ given analytically. We choose the 
form of the potential <& to be 

m = A=j, +1/2 (£), (20) 

with J £+1/2 being the semi- integer Bessel function of the first kind. In this case, 
the only modes in (|19|) correspond to £ = 2 and m = 2. The source n computed 
in this way determines the right-hand-side in Eqs. © and (fT0|) . The numerical 
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Test 1 




Radius \ Radius \ 

Figure 3. Comparison between the numerical and analytical solutions of Test 
1. Left figure shows the solution at the end of the run (r = 7.7) and right figure 
the absolute difference between the numerical and the analytical solutions in the 
inner region. 



solution obtained by the application of the present scheme is then compared to the 
corresponding analytical solution. The latter is expressed, in the inner domain, as a 
linear combination of the exact particular solution (|20p and the regular homogeneous 
solution of the operator (flT?]) . which can be expressed in terms of the semi-integer Bessel 
function of the first kind r\{£) — Ji+\/2 (£)) that can be easily evaluated because it 

reduces to the sum of £ + 1 terms. Regarding the homogenous solution in the external 
domain, this is expressed in terms of a second order Bessel function H£ +1 / 2 (£,)- Once 
the matching conditions at i?» and the Sommerfeld boundary conditions at £d are 
imposed, the explicit comparison with the analytical solution can be done, finding a 
very good agreement (cf. Figure [3] for a quantitative assessment). 

The tests above concern the case £ = 2 and m — 2, when a helical symmetry is 
present in the source. The number of points used are N v = Ng = 4, iV" mcr = 65 in 
the inner domain (for £ S [0,1]), 7V° utor = 257 in the outer domain (for £ e [1,21]), 
and a time step St = (1024) -1 . This case I = 2 and m = 2 has been chosen because 
it is relevant for gravitational wave emission in General Relativity. The choice of a 
symmetry in the source simplifies the comparison between analytical and numerical 
results. 

We have studied the convergence behavior of the code with the same settings, but 
changing the number of radial coefficients in both domains, while keeping fixed the 
ratio 7V° uter /iV™ ner = 8. In particular, we have checked that the difference between the 
numerical and analytical solutions would decay as e~ Nr , which is the expected rate for 
spectral approximation series of smooth functions. This check is displayed in Figured! 
where the discrepancy saturates at about 10~ 6 level, due to numerical limitations. 
There are two sources accounting for this accuracy level, whose behaviors were clearly 
identified in numerical investigations: first, the time finite-differencing errors, which 
are linked to the second-order time-scheme, and then the use of approximate boundary 
conditions (we remind that Sommerfeld boundary conditions are exact only for £ = 
modes) at £ = £d. We have checked that these latter errors would decay as £^ 2 , which 
is the predicted rate (see also [27]). 
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Test 1 



Absolute error in both domains 




40 



Number of coefficients (inner domain) 



Figure 4. Behavior of the accuracy obtained in Test 1 (see also Figure [3} as 
a function of the number of radial spectral coefficients. The displayed numbers 
correspond to the inner domain; in the outer domain eight times more coefficients 
have been used. 



4-2. Conservation of helical symmetry 



The second test to our resolution scheme studies how a solution $ to Eqs. |T]) and 
"tracks" a symmetry present in the prescribed matter source. More specifically, we 
consider a source n which is rigidly rotating with respect to the incrtial frame with 
an angular velocity f2. We choose its form as 

4 / \ 2" 



n 



'cos(2<^ + 2fit). 



(21) 



This source is symmetric with respect to the helical vector £ = dt — Qdip, i.e. 

Because of the existence of this symmetry, we expect the solution $ to present also 
the helical symmetry once a sufficient time has been elapsed (once it has reached some 
"stationary" regime) . 

We apply to this simple case the scheme described in section |3l Regarding the 
inner rotating domain, since the only non-vanishing terms correspond to i = 2 and 
m = 2, Eqs. ([5]) and (fTT))) read as follows: 



2„< 



d 2 a 



dr 2 
-C 2 



+ (2tt) 4 



db* 



&e + 

-(27T) 



dr 
2d_ 

da< 
1h 



(27r)4a 



n(0 



4(271-)^ 



2 d 



G 



(23) 



(24) 
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where for shortness, we have written respectively a < = a 22 , b < = b 22 and n = 
R 2 ■ C 2 n 22 . Note that in the inner rotating frame, the source n does not depends 
on r and h s = 0. In the exterior non-rotating domain we solve the similar expressions 
coming from Eqs. (jlip and (fT2")l . 

Initial conditions are chosen such that a = 6 = OatT = To, and the evolution 
is induced by the source n (see Fig. UK, where the profile at the first time-step is 
given). This was computed in a such way to ensure the continuity of the solution and 
its derivative across the boundary r — i?*. The dynamical evolution began at time 
r = t = 0, and the numerical settings were slightly different from those of the first 
tests, namely: A^ nner = 33, 7V° utor = 129, Ng = 4, N v = 4 and St = (512)' 1 . 

We study the emergence of the helical symmetry in the evolving dynamical 
solution in two different manners: 

i) First we assess the verification of the symmetry conditions given by the Eq. (|22|) . 
for the case of the constructed numerical solution. This provides a good self-consistent 
test of the accuracy of the solution. In the present case, we compute the quantity 5a, 
defined as 



where a is considered as a function of £ and r, and we study its evolution in time. 
Figure Oa shows the evaluation of (f2"5|) on the numerical solution: after the transient 
due to the switch on of the evolution, the solution must asymptotically verify da = 0. 
Relative error, computed at £ = 1 in the outer region, is plotted versus the time r. 
FigureEJb shows the difference between the computed solution and the solution of the 
Poisson equation for $ (instead of the wave equation |T]), to show propagation effects). 
The actual solution differs qualitatively with respect to the non-propagating one by 
a phase- (n/2) quadrature term that is showed in Figure [5]c. This term is responsible 
of the radiation reaction force (or braking force) , and we have checked that this term 
scales as (i?*/r) 5 as predicted by the analytical properties for an £ = 2 propagation 
mode. Note that the solution of the Poisson equation, i.e. satisfying the equation 
Aa < = — n, is 



where a = (1/20 - 1/14) and = (1/45 - 1/35). 

Finally, Figure [S] presents some radial profiles at various moments of the numerical 
solution to Eqs. and (|2"4")l. 

it) A second way to proceed in order to study the emergence of the helical symmetry 
in the solution, consists in comparing the previous dynamical solution with the one 
corresponding to a system obtained by neglecting the time derivative terms in Eqs. ([9])- 
(| 10(1 . in the inner domain. In this case the problem in the rotating region becomes 
an elliptic one, and it can be solved easily by inverting the operator L given in Eq. 
(|19|. whereas we still solve for the wave operator in the outer domain. Note that it 
is crucial for the ellipticity of this operator L the fact that the rotating grid is set 
inside the light cylinder. If this were not be the case, we would have a mixed-type 
differential operator (see e.g. [H]). Given the second-order Crank- Nicolson scheme we 
have employed, results computed with the method i) and the method it) differ only by 
quantities of second order in (St) 2 . As a complementary test to the emergence of the 
helical symmetry, Figure [7] shows evaluation of (|2"5"|) on the numerical solution to the 
symmetry- reduced operator (| 19|) . Note that the solution reaches the stationary regime 




(26) 
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Figure 5. Fig. [5^,: Relative error 8a for Test 2 in a logarithmic scale versus the 
time t, for Rl/R* = 3. The huge spike at time r = 1 is due to terms containing 
first time derivatives in Eqs. 1(9} , l|10|l that are switched on only at r = 1 in order 
to avoid a strong oscillation at the beginning of the simulation. The plateau 
reached by the error is due to time-differentiation and scales like St 2 . Fig. |5p: 
The actual solution of the wave equation in the rotating region (upper curve) and 
the solution of the Poisson equation. Fig. [5ji: Component of the solution <E> in 
quadrature with the source. 



much faster than when solving the full wave equation. The big bump appearing at 
the time r ~ 3.2 and the small one at the time r ~ 5.5 are generated by a partial 
reflexion at the external boundary of the external domain due to the approximate 
outgoing Sommerfeld boundary conditions [which are exact for i — 0; see Eq. ([13"])]. 
We have checked that the amplitude of the bumps scales as the inverse of the square 
of the dimension of the external grid, as predicted by the theory. 

5. Discussion and conclusions 

The rationale of this work is the following: a class of physical problems (rotating 
stars or orbiting stars/black holes) are easily treated in corotating grids. However, 
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Figure 6. Figures [6]a-c show some snapshots of the solution in Test 2, with 
Rl/R* = 2 at different times. At £ = 1 the "star" shows the frontier between 
the rotating and inertial domains. The vertical line at £ = 2 indicates the light 
cylinder. On Fig. |6]e the solution in the external region was multiplied by £ in 
order to better visualize the solution at the edge of the external grid. Figure [Blf 
is a zoom of the solution inside and near the light cylinder. The dimensionless 
"light velocity" is C = An. 
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Figure 7. Test of the helical symmetry 8a as in Fig. [5ji. 



analytical and numerical problems show up associated with the existence of a light 
cylinder. We have described a technique to match the solution computed in an inner 
rotating domain of external radius i?*, with the solution in an outer non- rotating 
domain under the condition i?„ < Rl, where Rl = cj Q is the light cylinder radius. 
The matching technique is independent of the manner of obtaining the solutions in each 
domain. In typical physical situations, the restriction R* < Rl on the radius is not 
severe: for example in the worst case for orbiting neutron stars, we have fl ~ 4000 s _1 
[29] and i?» ~ 20 km, providing a ratio Rl/R* ~ 4. As a more extreme case, this ratio 
for orbiting black holes before the merger is not smaller than 2. A typical problem to 
be treated with the technique described here, is the case two orbiting NS, for which 
the parameters of the orbit are known, and for which we want study small amplitude 
motions inside each star (convection r- modes and/or elliptical cylinder instabilities). 
The numerical study of these phenomena is easier in a domain in which the surface of 
the star does not move and boundary condition can be exactly imposed, as compared 
with a star traveling through the grid. 

Regarding the numerical implementation, this is not the only scheme making 
use of multi-domain spectral methods in order to deal with the resolution of partial 
differential equations in rotating domains (cf. |10j in the general dynamical case and 
[25] in the helically symmetric situation). The specificity of our approach is the use of 
a rotating domain (our inner domain) that, unlike the case in other schemes, does not 
extend up to infinity. In particular our inner rotating domain does not go beyond the 
light cylinder, thus providing a simple treatment of this analytical and numerical issue. 
The results we have obtained are accurate, and quantities like the back-reaction force 
are correctly computed. We conclude that our "stroboscopic" matching technique is 
a simple and efficient tool for the description of relativistic rotating systems, and in 
particular for avoiding analytical and numerical problems linked with the existence of 
a light cylinder. 

The technique described can be extended, mutata mutandis, to vectorial or 
tensorial equations, in particular choosing some convenient triad (or tetrad depending 
on the formulation), and projecting onto it the tensor components in order to deal 
with scalar quantities. 
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Finally, the case involving a non-constant angular velocity can also be treated: if 
the variation of the rotation period P in time is known, we can perform the coordinate 
transformation t = t(r) generalizing Eqs. ([3]), in such a way that wave equations 
correcting Eqs. ^ and ([3]) with extra terms in ^ are obtained. These additional 
terms are harmless for our scheme, provided that the ratio Rl / R* is sufficiently large. 
All those possibilities are yet to be explored. 
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